function [cost_measure spikeout mfr all_target_mfr] = spnet( conn_data_in, optional_pattern, all_target_mfr)

if nargin > 0
    if size(conn_data_in,2)==1
        conn_data_in = reshape(conn_data_in, 6, 3);
    end
end
start_save_sec = 0;
stop_sec = 1;


if nargin > 1
    r=optional_pattern;
else
     r=[.122 .33 .12 .47 1 .42 .03 .43 .31 .31 .45; ...
     0.2848    0.1522    0.1166    0.2534    0.4529    0.4443    0.1073    0.1756    0.2698    0.3289    0.7500; ...
     0.1543    0.2126    0.5000    0.1736    0.1401    0.0113    0.0370    0.2063    0.2159    0.2477    0.1655; ...
     0.1009    0.1317    0.1265    0.0578    0.0548    0.2500    0.1192    0.0408    0.0170    0.0298    0.1115];
%r=[    0.0321    0.4475    0.2013    0.9687    0.1278    0.2822    0.2924    0.1789    0.0484    0.2318    0.7734;...
%    0.7965    0.8375    0.8627    0.0655    0.6800    0.4220    0.9668    0.1903    0.9100    0.5607    0.0139;...
%    0.0341    0.7960    0.2558    0.2295    0.0283    0.6533    0.7794    0.0992    0.1220    0.5733    0.2672;...
%    0.0185    0.9369    0.0597    0.0331    0.8519    0.9216    0.0863    0.9099    0.3537    0.3028    0.9321];

    all_target_mfr = [...
         0.0007    0.0091   14.0048   37.5236   47.7772   37.5828   13.5727    0.0255    0.0067    0.0053    0.0084;...
         29.3920    7.6146    0.0011    0.0028    0.0073    0.0077    0.0003    0.0069   18.3303   35.9531   42.0335;...
         18.9614   25.8034   24.2515   11.4287    0.0012         0         0    0.0032    0.0051    0.0224    5.1803;...
           0.0009    0.0068    0.2328    4.4247    7.4911    8.7554    5.4898    0.3377         0         0    0.0006 ...
            ];


        
    epoch_ms = 1000;
    npatterns = 4;
    isi_ms = floor(epoch_ms/npatterns);

        
end

search_params = [1 3 4];
conn_data = zeros(6,9);
    

conn_arg = 1;

%function spikeout = spnet( stop_sec, start_save_sec,  conn_data)
% spikeout will be all spikes between start_savfe_sec and stop_sec.

global s sd syntype D N syndata v I Ne IN A AI vhist ehist ampahist nmdahist gabaahist gababhist Achand sum_tempi sum_tempe

draw = false;
learning = false;
  

spikeout=[];

% spnet.m: Spiking network with axonal conduction delays and STDP
% Created by Eugene M.Izhikevich.                February 3, 2004
% Modified to allow arbitrary delay distributions.  April 16,2008
% Rewritten by Jeffrey L. McKinstry, November, 2009 to allow anatomy.

D=5;                  % maximal conduction delay 
% excitatory neurons   % inhibitory neurons      % total number 
Ne=22;                Ni=22;                   N=Ne+Ni;

% set up group sizes.
IN = 1:11;
A = max(IN)+1:Ne;
Achand = max( A ) + 1:floor(Ne+Ni/2);
AI = max(Achand)+1:N;
N = max(AI);


% neuron parametes:
a=[0.02*ones(Ne,1);    0.1*ones(Ni,1)];
d=[   8*ones(Ne,1);    2*ones(Ni,1)];

nmda_ampa_ratio=ones(N,1);
nmda_ampa_ratio(Achand)=.2;
nmda_ampa_ratio(AI)=.12;

gabab_gabaa_ratio=0.0*ones(N,1);
gabab_gabaa_ratio(Achand)=0;
gabab_gabaa_ratio(AI)=0;

% Take special care not to have multiple connections between neurons
s = cell(D,1);
sd = cell(D,1);
syntype = cell(D,1);
for i = 1:D
    s{i} = sparse(N,N);
    syntype{i} = sparse(N,N);
    sd{i} = sparse(N,N);
end



% connections
%r=rand(1,12)*.5;  r(5)=1

patterns{1} = r(1,:);
patterns{2} = r(2,:);
patterns{3} = r(3,:);
patterns{4} = r(4,:); %[0 0 0 0 0 0 0 0 0 0 0 0];



%patterns{1} = [0 0 0 .5 1 .5 0 0 0 0 0 0];
%patterns{2} = [0 0 0 .5 1 .5 0 0 0 0 0 0];
%patterns{3} = [0 0 0 0 0 0 0 0 0 0 0 0];
%patterns{4} = [0 0 0 0 0 0 0 0 0 0 0 0];

%patterns{1} = [1 1 1 1 1 1 0 0 0 0 0 0];
%patterns{2} = [0 0 0 1 1 1 0 0 0 1 1 1];
%patterns{3} = [1 1 1 0 0 0 1 1 1 0 0 0];
%patterns{4} = [0 0 0 0 0 0 1 1 1 1 1 1];


syndata = [];


i=1;


radius=length(A)/2;
conn_data(i,:)=[0.5*radius 0 5 radius 3 3 .05 0 10]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(IN,A,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));


radius=length(AI)/2;
conn_data(i,:)=[0.5*radius 0 5 radius 3 3 .05 0 15];
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(IN,AI,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

radius=length(AI)/4;
conn_data(i,:)=[0.5*radius 0 2 radius 1 2 .05 0 15]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(A,AI,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));


radius=length(Achand)/4;
conn_data(i,:)=[radius 0 2 radius 1 2 .05 0 10]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(A,Achand,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

% make chandelier cells excitatory.
radius = length(A)/6;
conn_data(i,:)=[0.5*radius 0  1 radius 1 2 0.0 0 6];
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(Achand,A,'gauss', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));

radius = length(A)/6;
conn_data(i,:)=[0.5*radius 0 -6 radius 1 1 0.0 -6 0]; 
if nargin >= conn_arg,
    % merge the supplied parameters with defaults.
    conn_data(i,search_params)=conn_data_in(i,:);    
end
conn = conn_data(i,:);i=i+1; 
connect(AI,A,'surround', 'std', conn(1), 'initoffsetw', conn(2), 'initrangew', conn(3), 'radius', conn(4), 'delaymin', conn(5), 'delayrange', conn(6), 'lrate', conn(7), 'minw', conn(8), 'maxrangew', conn(9));


mfr=zeros(1,N);




%s=[6*ones(Ne,M);-5*ones(Ni,M)];         % synaptic weights

ltp_decay = exp(-1/20);
ltpcurve = .1*cumprod([ 1 ltp_decay*ones(1,49)]);
ltd_decay = exp(-1/20);
ltdcurve = -.12*cumprod([ 1 ltd_decay*ones(1,49)]);
stdpcurve = [ fliplr(ltpcurve)  ltdcurve 0];
  
if draw
    mv = zeros(N,1); % for visualization
    vhist = zeros(N,1000);
    ihist = zeros(N,1000);
    ehist = zeros(N,1000);
    ampahist = zeros(N,1000);
    nmdahist = zeros(N,1000);
    gabaahist = zeros(N,1000);
    gababhist = zeros(N,1000);
    itothist = zeros(N,1000);    
end


% restart for each simulation.
%rand('seed',1)
tau_ampa  = exp(-1/5);
tau_nmda  = exp(-1/150);
tau_gabaa = exp(-1/6);
tau_gabab = exp(-1/150);


v = -65*ones(N,1)+rand(N,1);                      % initial values
u = 0.2*ones(N,1);                             % initial values

MAXSPIKESPERSEC= N/2*1000;
firings=zeros(MAXSPIKESPERSEC,2);

firehist=cell(N,1);

I=zeros(N,1);
Igabaa = I;Igabab = I; Iampa = I;Inmda = I;
Ie=zeros(N,D);
Ii=zeros(N,D);

numspikes = 0;

snew = [];
for i=1:D
    snew = [snew s{i}];
end
snew=full(snew);
se=snew(1:max(Achand),:);
si=snew(max(Achand)+1:N,:);


for sec=1:stop_sec                      % simulation of 1 day
  for t=1:1000                          % simulation of 1 sec
    
    %randindex = ceil(N*rand);
    %I(randindex,1)=I(randindex,1)+20;                 % random thalamic input 
%    I(IN,1) = I(IN,1) + 20*patterns{ 1+floor((t-1)/250)};
    
    fired = find(v>=30);                % indices of fired neurons
    efired = v'>=30;
    ifired = efired(max(Achand)+1:N);
    efired = efired(1:max(Achand));

    v(fired)=-65;  
    u(fired)=u(fired)+d(fired);
    
    % deliver all the synapses to post targets with correct delays.
    % I is a shift register for each neuron.
    %for i = 1:D
     %  Ie(:,i) = Ie(:,i) + (spfired(1:max(Achand))*s{i}(1:max(Achand),:))';
     %  Ii(:,i) = Ii(:,i) + (spfired(max(Achand)+1:N)*s{i}(max(Achand)+1:N,:))';
    %end
     %Ie = Ie + reshape( (spfired(1:max(Achand))*snew(1:max(Achand),:)) , N,D);
     %Ii = Ii + reshape( (spfired(max(Achand)+1:N)*snew(max(Achand)+1:N,:)) ,N,D);
     Ie = Ie + reshape( efired*se , N,D);
     Ii = Ii + reshape( ifired*si ,N,D);


    Igabaa = Igabaa*tau_gabaa + Ii(:,1);
    Igabab = Igabab*tau_gabab + gabab_gabaa_ratio.*Ii(:,1);

    Iampa = Iampa*tau_ampa + Ie(:,1);
    Inmda = Inmda*tau_nmda + nmda_ampa_ratio.*Ie(:,1);

    xx = (-80-v)/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    
    I = Igabaa + Igabab + NMDAgate.*Inmda + Iampa;
    I(IN) = I(IN) + 40*patterns{ 1+floor((t-1)/isi_ms)}';
        
    v=v+0.5*((0.04*v+5).*v+140-u+I(:,1));    % for numerical 
    v=v+0.5*((0.04*v+5).*v+140-u+I(:,1));    % stability time 
    u=u+a.*(0.2*v-u);                   % step is 0.5 ms
    
    if draw, ihist(:,t) = Ii(:,1); ehist(:,t)= Ie(:,1); 
        itothist(:,t) = I;
        ampahist(:,t)= Iampa; nmdahist(:,t)= NMDAgate.*Inmda; 
        gabaahist(:,t)=Igabaa; gababhist(:,t)=Igabab;
    end
    
    %I= [I(:,2:D) zeros(N,1)];  % rotate.     
    Ie= [Ie(:,2:D) zeros(N,1)];  % rotate.     
    Ii= [Ii(:,2:D) zeros(N,1)];  % rotate.     
    
    
    % save the spike data.
    firings(numspikes+1:numspikes+length(fired),:)=[t*ones(length(fired),1),fired];
    numspikes = numspikes + length(fired);
    
    if learning
        for i =1:length(fired)
            firehist{fired(i)}(end+1)=t;
        end
    end
    %imagesc(I);
  
    if draw
        vhist(:,t)=min(30,v);
        
        tau = exp(-1/200);
        mfr(fired)=mfr(fired) + (1-tau)*1000;
        mfr=mfr*tau;

        mv = tau*mv + (1-tau)*v;
        if mod(t,250)==0
            figure(3);
            bar(mfr)
            drawnow;
            
            %figure(4)
            %plot(mv);
            %drawnow;
        end
    end
    
  end;
  
  
  disp(['t=' num2str(sec)])
  
  if draw
      
      figure(1)
      plot(firings(:,1),firings(:,2),'.');
      axis([0 1000 0 N]); drawnow;

      figure(2)
      hist(s{1}(find(s{1})),50);

      figure(4)
      subplot(7,1,1)
      imagesc(vhist,[-70 -60]); colorbar;
      title('vhist')
      subplot(7,1,2)
      imagesc(ihist); colorbar;
      title('inhib input history');
      subplot(7,1,3)
      imagesc(ehist); colorbar;
      title('excit input history');
      subplot(7,1,4)
      imagesc(ampahist); colorbar;
      title('ampa input history');
      subplot(7,1,5)
      imagesc(nmdahist); colorbar;
      title('nmda input history');
      subplot(7,1,6)
      imagesc(gabaahist); colorbar;
      title('gabaa input history');
      subplot(7,1,7)
      imagesc(gababhist); colorbar;
      title('gabab input history');
      
      figure(5)
      subplot(2,1,1)
      j=1;
      for i = 1:100:599
          sum_temp(:,j) = sum(itothist(:,i:i+99),2)  ;
          j=j+1;
      end
      plot(sum_temp);
      size(sum_temp)
      
      legend('1','2','3','4','5','6');
      subplot(2,1,2)
      j=1;
      for i = 1:100:599
          sum_tempi(:,j) = sum(ihist(:,i:i+99),2) ;
          sum_tempe(:,j) = sum(ehist(:,i:i+99),2) ;
          j=j+1;
      end
      plot(sum_tempi);
      size(sum_temp)
      
      
  end
  % do some learning based on spike times.
  % sd in Eugene's spnet did not get used until after 1 second, so this is
  % equivalent.  Just makes for a cleaner implementation.
  
  
  if learning
      for delay = 1:D
        % get all the nonzero elements from synaptic matrix.
        [pre_i post_i] = find( s{delay} );

        % Adjust these weights.
        for j = 1:length(pre_i)
            type = syntype{delay}(pre_i(j),post_i(j));
            w=s{delay}(pre_i(j),post_i(j));

            %if type == 6, w,  sd{delay}(pre_i(j),post_i(j)),end

            if syndata(type).max < .00001
                scale = 1;
            else
                scale = w/ syndata(type).max;
            end
            sd{delay}(pre_i(j),post_i(j)) = sd{delay}(pre_i(j),post_i(j)) + stdp(firehist{pre_i(j)},firehist{post_i(j)}, delay, stdpcurve,syndata(type).lrate, scale);

            %if type == 6, w,  sd{delay}(pre_i(j),post_i(j)),end
            s{delay}(pre_i(j),post_i(j)) = max( syndata(type).min ,min(syndata(type).max, syndata(type).lrate*.01 + w + sd{delay}(pre_i(j),post_i(j)))); 
        end

        sd{delay}=0.9*sd{delay};

      end
      % keep spikes near the seconds boundary for future STDP calculations.
      for i = 1:N
        firehist{i}=firehist{i}( find(firehist{i}> 950) )-1000;
      end
  end
  
  if sec >= start_save_sec
      spikeout = [spikeout; [firings(1:numspikes,1)+1000*(sec-1) firings(1:numspikes,2)]];
  end
  
  numspikes = 0;
  
  
end;

% Calculate cost measure

% build the desired mfr profile.


% all_target_mfr = zeros(npatterns,length(A));
% for pat = 1:npatterns
% 
%     std=length(A)/10;
%     radius=floor(length(A)/2);
%     target_mfr=(gaussian_1d(std,-radius:radius));
%     target_mfr=target_mfr*(1/max(target_mfr)).*30; % winner fires at 30 hz.
%     [m max_index]=max( patterns{pat} );  % shift s.t. winner is in the middle.
%     target_mfr=target_mfr(1:length(A));
%     % assumes length(A) = length(r);
%     mid = ceil(length(A)/2);
%     shift = abs(max_index - mid)
%     if max_index > mid
%         % shift right
%         target_mfr=[ target_mfr(end-shift+1:end) target_mfr(1:end-shift)];
%     else
%         % mid >= max_index
%         target_mfr=[ target_mfr(shift+1:end) target_mfr(1:shift)];    
%     end
% 
%     all_target_mfr(pat,:)=target_mfr;
% end
    
% bin the spike data by pattern.
if ~isempty(spikeout)    
    S=sparse( spikeout(:,2), spikeout(:,1), 1, N , stop_sec*1000 );

    % Gather mean rates for stats.

    t=start_save_sec*1000 + 1;
    count = zeros(1,npatterns);
    mfr = zeros(npatterns,N);
    for i = start_save_sec*1000:epoch_ms:stop_sec*1000-1
        for p = 1:npatterns
            mfr(p, : ) = mfr(p, :) + (1000/isi_ms)*full(sum( S( :, t+20:(t+isi_ms-1)), 2))';
            count(p)=count(p)+1;
            t=t+isi_ms;
        end
    end
    mfr=mfr./(repmat( count', 1, size(mfr,2) ) + .000001);


    % Look at sparsness measures
    %     
    %     j= 1;
    %     for i = max(IN)+1:Ne
    %         lifetime_sparseness(j) = sparseness( mfr(:,i), 40 ); % 40 Hz is max
    %         j= j + 1;
    %     end
    %     disp(['Mean lifetime sparseness: ' num2str( mean(lifetime_sparseness))]);
    %     
    %     for i = 1:patterns
    %         population_sparseness(i) = sparseness(mfr(i,max(IN)+1:Ne), 40);
    %     end
    %     disp(['Mean population sparseness: ' num2str( mean(population_sparseness))]);

        % save the best
    %    d= sqrt( (mean(population_sparseness)-0.5)^2 )
    %d= sqrt( (mean(population_sparseness)-0.5)^2 + (mean(lifetime_sparseness)-0.5)^2)

    cost_measure = sqrt( sum(sum((mfr(:,max(IN)+1:Ne) - all_target_mfr ).^2 )))

else
    cost_measure = 1e10;
end





%
% possible arguments:
%  arbor_type == 'gauss'
%  'std', 'minw', 'maxrangew', 'dist', 'delayrange', 'lrate';
%
function connect( from_indices, to_indices, arbor_type, varargin)

global s sd syntype N D syndata


L = length( syndata);


for k= 1 :2: size(varargin,2) 

        switch lower( varargin{k} )
          case {'std'}
            std = varargin{k+1};
          case 'initoffsetw'
            minw = varargin{k+1};
          case 'initrangew'
            maxw = varargin{k+1};
          case 'radius'
            radius = round(varargin{k+1});
          case 'delaymin'
            delaymin = round(varargin{k+1});
          case 'delayrange'
            delayrange = round(varargin{k+1});
            if delaymin + delayrange-1 > D
                disp('delay > max delay in connect. Using max.');
                delayrange
                delaymin
                delayrange = D-delaymin + 1;
            end
          case 'lrate'
            syndata(L+1).lrate = varargin{k+1};
          case 'minw'
            syndata(L+1).min = varargin{k+1};
          case 'maxrangew'
            syndata(L+1).max = syndata(L+1).min + varargin{k+1};
          otherwise
            disp(['Unknown param.' varargin{k}])
        end

end

% Assume wraparound.

% clip radius to not touch target neurons more than once.
if radius > length(to_indices)/2-.5, radius = floor(length(to_indices)/2-.5);end

p=[to_indices to_indices to_indices]; %randperm(N);
% create the list of indices for each neuron; project topographically.
step = length(to_indices)/length(from_indices);
rel_index = .5;

if strcmp(arbor_type, 'gauss')
    w = minw+maxw *std*sqrt(2*pi)* gaussian_1d(std, -round(radius):round(radius) );
else 
    % surround
    x=-round(radius):round(radius);
    w=std*sqrt(2*pi)* abs(x).*gaussian_1d(std, x );
    w = minw+maxw*(1/max(w))*w
end
for i=from_indices
    my_center = round(rel_index*step);
    rel_index=rel_index + 1;
    myoffset = length(to_indices) + my_center;

    target_indices = p( myoffset + (-round(radius):round(radius)) );
    
    delays = delaymin + floor( delayrange*(abs(-round(radius):round(radius))/(round(radius)+.001)));

    for j = 1:length(target_indices)
        s{delays(j)}(i,target_indices(j)) = w(j); % + noise.
        syntype{delays(j)}(i,target_indices(j)) = length( syndata );
        sd{delays(j)}(i,target_indices(j)) = 0.0000001; % + noise.        
    end

    
end;





